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the piston to the expansion chamber. It is then possible to simulate the dynamic (steady state stroke 
and operation frequency) as well as the thermodynamic performances (output power and efficiency) 
for given mean pressure, heat source and heat sink temperatures. The motion amplitude especially can 
be determined by the spring-mass properties of the moving parts and the main nonlinear effects which 
are taken into account in the model. The thermodynamic features of the model have then been validated 
using the classical isothermal Schmidt analysis for a given stroke. A three-phase low temperature differ¬ 
ential double acting free membrane architecture has been built and tested. The experimental results are 
compared with the model and a satisfactory agreement is obtained. The stroke and operating frequency 
are predicted with less than 2% error whereas the output power discrepancy is of about 30%. Finally, some 
optimization routes are suggested to improve the design and maximize the performances aiming at 
waste heat recovery applications. 
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1. Introduction 

Low temperature differential (LTD) Stirling engines are ma¬ 
chines that can operate with hot source temperature of about 
150 °C. At this operating temperature and assuming a sink temper¬ 
ature of 25 °C, the Carnot efficiency is 29%. Then, a Stirling gener¬ 
ator which would achieve 50% of this efficiency might be of 
interest for power generation. Solar powered LTD Stirling engine 
appears to be a promising technology 1 ] and is also a potential 
technology for waste heat recovery [2,3], In both applications, 
the simplicity and reliability of the Stirling machines are significant 
advantages toward the development of low cost generators. 

Among the various potential Stirling architectures, the double 
acting engines type also called Rinia or multiphase architecture 
[4,5] is of particular interest for the aforementioned applications. 
In such configuration, three, four or more alpha type engines (also 
called phases) are connected to each other. The piston of the 
engine i-1 acts as the displacer of its neighboring engine i. In [6], 
Abdullah studied such a double acting LTD Stirling. He underlines 
that the specific material and components required for the 
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crankshaft, the connecting rods and the gaskets are an important 
issue for LTD applications. 

A free piston Stirling engine (FPSE) architecture allows avoiding 
complex mechanical linkages and the resulting robustness and 
reliability questions. Though, their optimization has been proven 
to be difficult especially for the piston-displacer phase angle con¬ 
trol [7,8], The double acting architecture overcomes this difficulty 
because the phase angle is fixed by the number of phases. It equals 
120° in the case of three engines and 90° in the case of four 
engines. Therefore, the free piston double acting arrangement 
appears to be a great advantage compared to the usual FPSE. 
However, a proper sealing of the piston and the displacer can be 
difficult to ensure. 

This last obstacle can be classically solved using membranes. In 
his work, Minassians proposed a LTD double acting FPSE using 
membranes instead of pistons [9,10]. The operation has been dem¬ 
onstrated on a 3 x 350 cm 3 total volume engine. A theoretical effi¬ 
ciency of 19.7% which is about 70% of the Carnot efficiency was 
expected to be reached using helium or pressurized air engine. 

In the case of double acting FPSE, the strong coupling between 
dynamic and thermodynamic has to be addressed properly. The 
major non-linear effects have to be integrated in the model to pre¬ 
dict the performances. These effects are usually associated to the 
gas friction losses within the components and the mechanical 
non-linearity [7,8,11-13], 
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Nomenclature 


A d displacer area (m 2 ) 

A v piston area (m 2 ) 

Aff free flow area (m 2 ) 

A w wetted area (m 2 ) 

c d displacer mechanical dissipation (NnT 1 s) 

c p piston mechanical dissipation (Nm -1 s) 

Cf friction factor (-) 

D damping coefficient (Nm -1 s) 

d w wire diameter of the matrix regenerator (m) 
d hi hydraulic diameter of component i (m) 

D d displacer diameter (m) 

D p piston diameter (m) 

E b young modulus (NnT 2 ) 

/ frequency (Hz) 

ff friction factor (-) 

f m mechanical force (N) 

I b second moment of the spring beam (m 4 ) 

j complex number such as j 2 =-1 

k gas thermal conductivity (WnT 1 I<) 

K spring beam stiffness (NnT 1 ) 

Lj length of exchanger (m) 

L b length of the spring beam (m) 

m tot total mass of gas (kg) 

M mass of a rigid link (kg) 

m, mass of gas within component i (kg) 

n, - number of channels of component i (-) 

P m mechanical power (W) 

p pressure (N nr 2 ) 

r ideal mass gas constant (J kg -1 K _1 ) 

S b cross-section of the spring beam (m 2 ) 

T cold cold temperature (K) 


T bot hot temperature (I<) 

T, temperature of gas within component i (K) 

T w wall temperature (K) 

u, mean gas velocity in component i (ms -1 ) 

Vi volume of component i (m 3 ) 

V sw swept volume (m 3 ) 

Vswe expansion space swept volume (m 3 ) 

Vswc compression space swept volume (m 3 ) 

Greek symbols 

a swept volume phase angle (rad) 

X fluid flow resistance coefficient (s -1 ) 

Ap pressure loss (N m -2 ) 

s scale (-) 

k swept volume ratio (-) 

/( dynamic viscosity (N nT 2 s) 

p gas density (kg nT 3 ) 

p b material density (kg nT 3 ) 

co pulsation (rad s -1 ) 

If porosity of the regenerator (-) 

Subscripts 
b beam 

c compression space 

d displacer 

e expansion space 

h heater 

k cooler 

p piston 

r regenerator 


In [9] Minassians proposed an isothermal dynamic model to 
study the behavior of the multiphase engine. Based on the Schmidt 
analysis, the instantaneous pressure is expressed as a function of 
the piston positions. Then the mechanical equilibrium equation 
of each piston is established including the dissipation forces and 
the external load. A linearization strategy allows an analytical 
expression of the start-up condition to be obtained. Above this crit¬ 
ical value, a numerical resolution has to be used to obtain the peri¬ 
odic motion of the engine and to evaluate its performance. The 
flow friction and hysteresis losses have been pointed out because 
of their effects on the engine performances and have to be taken 
into account for design of such LTD engines. 

Another modeling approach based on the control system theory 
of a wobble-yoke double acting Stirling engine has been presented 
in [14,15]. The engine architecture implies some flexibility in the 
mechanical connections of the pistons so its operation is close to 
double acting FPSE. The modeling strategy relies on the dynamical 
equilibrium equations established for each of the piston. The pres¬ 
sure can then be defined as a function of the pistons’ position using 
the Schmidt analysis [4] (perfect gas law, isothermal transforma¬ 
tions). A linearization is performed to obtain a state space repre¬ 
sentation of the system. Such an approach leads to an unstable 
dynamics so a control concept based modification namely a pre¬ 
compensator, is defined to reproduce the pistons’ amplitudes. This 
model has not been compared to experimental results and some 
important losses such as flow friction are not taken into account. 

In this paper, a generic network model is proposed and vali¬ 
dated to study the dynamic and thermodynamics behavior of a 
double acting free membrane Stirling engine. 

As a basis for the study, a novel configuration of a LTD double 
acting FPSE derived from the work of Minassians [10] is 


introduced. Fig. la presents the global engine. It is composed of 
three alpha type engines using membranes and beam springs 
which are connected to each other via rigid links. The top side of 
the system can be connected to the heat source whereas the bot¬ 
tom side is cooled. For each engine, the simplest architecture has 
been sought here. So, the heater and cooler, the expansion and 
compression chambers are identical (see Fig. lb). The heater and 
cooler consist in 20 holed punched plates fitted together in a total 
thickness L b = L k . By doing this, the axial conduction through the 
solid material is impeded whereas the radial conduction allows 
aiming at isothermal processes in the heat exchanger volumes. 
The regenerator consists in a stack of 40 plain weave wire meshes. 
The wire diameter is 50 pm and apertures are 350 pm squares. The 
main dimensions and characteristics of the engine are given in 
Table 1. 

In a first part, the equivalent electrical model for each compo¬ 
nent of the engine is detailed. The classical electrical analogy is 
chosen: the pressure is seen as the electrical voltage and the vol¬ 
ume flow rate as the electrical current. The continuity and momen¬ 
tum equations are established for the compression and expansion 
chambers, the heat exchangers and the regenerator. Then using lin¬ 
ear approximation for the time and space variables, the equivalent 
electrical circuits are derived. Finally, the electrical circuits are 
connected to each other to obtain the dynamic equations. A global 
analysis is then achieved to establish the start-up condition and 
evaluate the performances. Then, as a first validation step, the ther¬ 
modynamic results of the proposed model are compared with 
those obtained using the Schmidt analysis. The effect of the load 
on the critical start-up temperature and the engine performances 
is studied. An experimental prototype has been built and is 
presented in the last part. The experimental measurements are 
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Phase 1 


Phase 2 



Fig. 1. (a ) Schematic of the proposed double acting free-piston engine and (b) exploded view of a single phase. 


Table 1 

Main parameters values of the studied double acting Stirling engine. 


Volumes for one engine 

Vc 

164.4 cm 3 


164.4 cm 3 

Vtct 

440 cm 3 

Heat exchangers and regenerator 
Hydraulic 0 

4 

4 mm 

4 

4 mm 

dr 

0.28 mm 

Length 

4 

15.5 mm 

4 

15.5 mm 

L r 

15.2 mm 

Beam springs 

Material = Steel 

4 

100 mm 

S b 

5.6 mm 2 

i b 

0.146 mm 4 

Membrane 

Material = silicone 

A p 

4.3 cm 2 

Ad 

4.3 cm 2 



Typical temperatures and mean pressure 







T c 

25-40 °C 

T„ 

140-160°C 

p 

Atmospheric pressure 


compared with the proposed model and the discrepancies are 
studied to improve the prototype and appraise the modeling strat¬ 
egy. Based on the presented prototype, technological develop¬ 
ments are suggested to improve the design and enhance the 
performances aiming at waste heat recovery applications. 

2. Network modeling of each component 

In the following the components of one engine (or phase) are 
addressed successively. The governing equation are set and linear¬ 
ized to obtain an equivalent electrical model. 

2 . 1 . Compression and expansion spaces 


fnj( t )=^-Vj(t)+^-Pj(t) ( 2 ) 

The mass flow rate out of the considered volume is: 
rhj(t) = -pjUj(t) in which the local density is pj = T and tij(t) is 
the volumetric flow rate out of the chamber. u p (t) = -Vj(t) is the 
instantaneous volume swept by the membrane and is then positive 
whenever the volume decreases (Vj(t) < 0). Therefore, Eq. (2) is fi¬ 
nally expressed as: 

Uj(t) = u p (t) - Cjjjj(t) ( 3 ) 

In which Cj = ^ is an equivalent electrical capacitance. Eq. (3) 
can then be read as the electrical circuit of Fig. 2a. 


The compression and expansion spaces are the volumes swept 
by the membranes. An isothermal assumption is chosen here. 

The continuity equation for a perfect gas within a chamber 
leads to: 


2.2. Heat exchanger 

We assume that the inertia of the gas is neglected and the tem¬ 
perature is constant throughout the heat exchanger. The local 


rrtj(t) 


Pj(t) 
rTj(t ) 


Vj(t) + 


Vj(t) 

rTj(t) 


Pj(t) 


( 1 ) 


In which j = c, e for the compression and expansion spaces 
respectively. 

The variations of the pressure, the mass and the volume are as¬ 
sumed to be small with respect to their mean value over a cycle 
and are denoted p, fhj and Vj respectively. Then, we can set: 
Pj(0 =P+Pj(t), mj(t) = fhj + ifij(t), Vj(t) = Vj+Vj(t). Neglecting 
high order terms, Eq. (1) is written again as: 


u p u. 



(a) 


U 0 Ul Uq U Lr 



Fig. 2. Equivalent electrical circuits for: (a) the expansion or compression chamber, 
(b) the heat exchanger and (c) the regenerator. 
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equations for the conservation of mass and momentum for ID flow 


are then: 


1 dp q (x, t) 1 dm q {x, t)_ 
rT q dt A ffq dx 

(4) 


(5) 


In which q = k, h for the cooler and the heater respectively. 
Applying the same approximation as before and retaining the 


lowest order terms, we obtain: 


1 dp q (x, t) 1 dm q (x. t) 
rT q dt A ffq dx 

(6) 


(7) 


In the Laplace space, the previous Eqs. (6) and (7) can be inte¬ 
grated for the geometric variable x to bring: 


m q (x,s) 


cosh(T(s)x) 

-sW sinh < r (sW 

m,(0,s) 

_Pq(X,S) 


-Z c (s) sinh(T(s)x) 

cosh(r(s)x) 

Pq(0,S) 


( 8 ) 

In which T(s) = Z c (s) = ^\j\xrT q and s is the Laplace 

transformation variable. 


Retaining the lowest order terms of the power series expansions 
of the cosh and the sinh functions, a spatially independent model is 
derived from Eq. (8). The governing equations for the heat 
exchangers are then: 

m q (L q ,s) = rh q (0,s)-^kksp q (0,s) (9) 

11 q 

P q (L q ,s) =p,(0,s)-|^m,(0,s) (10) 

A ffq 

Replacing Eq. (9) in Eq. (10), the relationship between the input 
and output pressures can then be expressed as: 

p q (L q ,s) = ^ -s^fjp q (0,s) ~^m q (L q ,s) (11) 

Retaining the first order terms, Eq. (11) is finally simplified to: 

p q (L q ,s) =p„(0,s) -^m q (L q ,s) (12) 

A ffq 

Using the volumetric flow rate and the time variable, Eqs. (9) 
and (12) are rewritten as: 

Uu,(t) = u 0q (t)-C q p 0q (t) (13) 

PLq(t) = Poq(t) ~ ^PqU Lq (t ) (14) 

A ffq 

In which C q = is the equivalent capacitance and p q = ^ 

2.2.1. Modeling of the pressure drop 

The fluid flow resistance coefficient / can be derived from the 
pressure drop - Reynolds number relationship and is related to 
the friction factor/)by Eq. (15): 

Z = 8^ r |u t ,(t)| (15) 

In which/)is in the form/)= CfRey" in which Rey = 4^-^^(t) is 
the instantaneous Reynolds number. 


Using Eq. (15) in Eq. (14), we obtain the expression of the pres¬ 
sure drop across the heat exchanger: 

n2n+5 jjn ■ ! , 

Pu,(t) =p 0q (t) (16) 

Eq. (16) is nonlinear and it is useful to anticipate the value of the 
friction factor C f for the operating conditions of the engine for fur¬ 
ther simplification. 

LTD Stirling engines usually presents low operating frequency 
and small swept volume. Assuming millimeter range strokes as 
well as operating frequency below 100 Hz, the mean Reynolds 
number can be approximated using the model proposed by Organ 
[16 and used in [7], Fig. 3 presents the estimation of the Reynolds 
numbers in the heat exchangers and the regenerator for membrane 
strokes of 3 mm, T e = 130 °C, T c = 50 °C and a phase angle of 120°. 

Following Martini’s approximation for circular tubes [17], Rey¬ 
nolds numbers less than 2000 lead to C/= 16 and n = — 1. 

In these conditions, Eq. (16) can be simplified to: 

PLq{t) =Poq(t) ~R q U Lq (t) (17) 

In which R q = 128^- is the equivalent constant resistance. 

Finally, Eqs. (13) and (17) can then be seen as the electrical cir¬ 
cuit of Fig. 2b. The evaluation of the Reynolds number from the 
model results will be performed and this assumption validated. 


2.3. Regenerator 


Following the same strategy as for the heat exchangers, the gov¬ 
erning equations for the regenerator are: 


m r (L r ,s) = m r (0,s) - J l ffrL ' sp r (0.s) 

• 1 rmean 


(18) 


Pr(L r ,S) =p r (0,S) ~^fh r (L r ,S) (19) 

A ffr 

In which T lrnean = i ^fr c r 

An electrical equivalent circuit is finally completed using the 
volumetric flow rate: 

U Lr (s)=^Uor-C r SPr(0,S) (20) 

1 C 


PLr(s) = Po(S ) - ^PrULr(s) (21) 

A ffr 

In which C r = f ^ is the equivalent capacitance. Eq. (20) 

can be written again considering an additional “current source” 
to take into account the temperature gradient in the regenerator: 

Ms) = Ms) - C r sp r (0, s ) + ( y - 1 ) Ms) (22) 



freq [Hz] 


Fig. 3. Estimation of the Reynolds numbers for the heat exchangers and the 
regenerator as a function of the operating frequency. 
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The last term of Eq. (22) is a thermal source flow which will be 
denoted u th (s) in the following defining u tb (s) = Cu 0r (s) and 
G = (pi - 1). Fig. 2c presents the equivalent electrical circuit. 

2.3.1. Modeling of the pressure drop 

As for the heat exchangers, the resistance coefficient % is related 
to the friction factor and the Reynolds number. For a woven screen 
matrix regenerator the hydraulic diameter is a function of the wire 

diameter and the porosity such as: d h = d w For the regenerator, 

the estimated Reynolds number is smaller than 30 as can be seen 
in Fig. 3. Then, the proposed experimental correlation is obtained 
from [18]: Cf =3.63 and n = -0.84. The dissipation term in Eq. 
(21) is finally a nonlinear expression with respect to the volumetric 
flow rate u Lr (t). 

A quadratic polynomial approximation is used so the equivalent 
resistance R r is defined as: 

1 t jji+l I 

Rr = 2 C /%TTOS (a(n) + b(n) Mt )2) (23) 

" “hr ™ffr 

In which a(n) and b(n ) are functions of the coefficient n. 

A numerical evaluation in the case of the proposed Stirling 
show that the quadratic term b(n) is several orders of magnitude 
lower than the constant term a(n). Therefore, only the latter will 
be retained in the simulation leading to a linear damping effect. 
Eqs. (21) and (22) can then be presented as the equivalent electri¬ 
cal circuit of Fig. 2c. 

2.4. Beam spring stiffness and rigid arm dynamics 

The dynamics of the engine partly relies on the beam spring and 
the mass of the rigid links which connect the membranes together 
(see Fig. la), whereas the silicone membranes mostly induce 
mechanical dissipations. 

2.4.1. Beam spring stiffness 

Beam springs are clamped to the rigid link at one end and to the 
frame of the engine at the other side. The stroke of the engine is 
bounded by the strongest nonlinearity of the engine [13], Then, 
large displacements of the beam are to be considered and may 
be the highest nonlinear effect because the pressure drop effects 
are expected to be small (low Reynolds numbers). 

Following the modeling strategy of Azrar [19], a quasi-static 
model based on the Lagrange’s equations is performed to obtain 
the large deformation stiffness of each spring beam. Using a single 
mechanical eigenmode approximation, the displacement of the tip 
of the beam can be written as: 

w(L b ,t) = q t (t)</>j(L b ) (24) 

The static equilibrium equation of one beam spring is given by 
Eq. (25): 

fc„g 1 (t) + 2b„„q 1 (t) 3 = F 1 / m (t) (25) 

In which, k u =E b I b f L 0 b ^^dx, b un =^f^^dx 
jL b d^jx) d^jx)^ = (p^L b ) and f m {t ) any additional force applied 
to the tip of the beam. 

From Eqs. (24) and (25), the nonlinear stiffness of the beam 
spring is defined as: 

I<(w(L b , t)) = k u + w(L b . tf (26) 

<M4) 

In order to validate the proposed semi-analytical approach, a 
FEM model of the beam has been performed and the stiffness 
evaluated. Fig. 4 shows the evolution of the nonlinear stiffness as a 
function of the tip of the beam displacement for the FEM and the 



Fig. 4. Evolution of the nonlinear spring beam stiffness as a function of the 
deflection for the beam dimensions given in Table 1. 

semi-analytical model. Though a moderate discrepancy can be seen 
between the two models, the semi-analytical model is validated and 
precise enough for the proposed engine modeling strategy. 

2.4.2. Rigid link dynamics 

The dynamic equilibrium of each rigid arm linked to the piston 
and displacer membranes of two successive phases can be written 
as in Eq. (27): 

M dw(L b ,t) + D dw(L b ,t ) + 2 , <(w(Lb t))w{Lb t) = ApPj j(f) _ Adp . {t) 
at ai 

(27) 

D is the damping coefficient which represents the mechanical losses 
within the membranes and any additional load (D = D mec + Di oad ). 

Recalling Eq. (3), the volumetric flow rate related to the dis¬ 
placement of the tip of the beam is u p (t) = - A p w(L b , t) for the pis¬ 
ton membrane (u d (t) = A d w(L b .t) for the displacer membrane). 
Moreover, because the piston and displacer areas are the same, 
Eq. (27) can finally be written as: 

LpdU 3r+ R p u p(t)+c p {u p ) / u pW =Pi(Q — Pi-i (0 ( 28 ) 

A 2 

In which L p = R p = ^ and C p (u p ) = are the equivalent 
inductance, resistance and nonlinear capacitance respectively. 

The electrical circuit of Fig. 5 is then consistent with Eq. (28). 
The nonlinearity of the stiffness term in Eq. (27) (i.e. the capac¬ 
itive term in Eq. (28)) as to be taken into account properly for the 
system analysis in the Laplace space. Numerous methods can be 
used to obtain analytical solutions for the dynamic of nonlinear 
oscillators as described in [20]. The approximate first harmonic 
solution of free oscillations of Eq. (28) when p,_i(t) = p,(t) = 0 can 
be written as: 

u p (t) = u p0 cos(co p t) = u „ O cos ^x/o7(l + | 0 u p O ) f j (29) 

In which oq = a 3 = 4 tlnl , -4 and u„ 0 is the initial ampli- 

™ (p,(L b ) MA p 

tude condition of the differential equation. 

Eq. (29) exhibits the characteristic frequency amplitude depen¬ 
dent behavior of nonlinear cubic oscillator. The natural oscillation 
frequency is a function of the amplitude of the response. Therefore, 



Fig. 5. Electrical circuit equivalent to the dynamics of the membrane and beam 
spring. 
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Fig. 6. Global network model of the engine i. 


we approximate the Laplace transform as written in Eq. (30) in In the Laplace space, the input-output relationships are then 

which the “natural” frequency is dependent on the amplitude of expressed in Eqs. (33) and (34): 
the volumetric flow rate U p : 


e [1 + sZ p (Cc + Cr) + G(1 + sCni?r)(l + sCcZp) + sCn(Rr + Zp + sR r Z p (Cc + C r ))]u,_i — u, 
' i “ 1 ” s[C H + C r + sC H C r R r + Cc( 1 + G)(l + sC H R r )] 


(33) 


(l p s 2 +R p s + 1^,(1 +||u p 2 ) 2 ju p = s(p,(s) -p,_,(s)) (30) 

2.5. Hysteresis loss 

The compression and expansion spaces act as gas springs and 
hysteresis loss is considered within these spaces as it was pointed 
out in Minassian’ study [9], Following the approximation of Urieli 
[4], the average value of the gas spring hysteresis loss power for 
the compression or the expansion space is: 

W lm = ^cof(y-l)T w pk(^j A w (31) 

In which AVj is the average swept volume. 

Therefore, a linear damping can be inferred from Eq. (31) and is 
defined as [4]: 

D hyst = \l~yi(y--l)T w pk^A w (32) 

This damping effect can be added to the mechanical dissipation. 
Consequently, the equivalent resistance R p of Eq. (28) is modified 
such as the hysteresis loss is taken into account and the equivalent 
electrical resistance is: R p = 4 + D hyst . 

A p 

3. Global model: coupled thermodynamic-dynamic analysis 

The modeling of each of the elements of one free piston engine i 
has been performed. Using the continuity of the pressure and the 
volumetric flow rate from one element to the other, a global elec¬ 
trical network is derived and is shown in Fig. 6. The pressure pf_, at 
the left hand side of the scheme stands for the expansion chamber 
pressure of the engine i - 1 whereas p? is the expansion chamber 
pressure of the considered engine i. 

From this equivalent network, the relationships between the in¬ 
put pressure pf_, and the volumetric flow rate u,_! and the output 
pf and Uj can be established. 

In the following, the pressure drops through the heat exchang¬ 
ers is neglected compared to the one across the regenerator. There¬ 
fore, the equivalent electrical resistances R k and R h are supposed to 
be zero and the scheme of Fig. 6 is simplified accordingly and is 
presented in Fig. 7. 


(1 + G)u i _ 1 - (1 + s(C r R r + CcQ + G)R r ))Uj 
Pi s[C„ + C r + sC H CA + C c (l+G)(l+sC H J? r )] 1 J 

In which Z p = ^ + R p + sL p , C c = C c + C k and C H = C e + C h . 

3A. Solution procedure 

For each of the engine, two equations similar to Eqs. (33) and 
(34) are obtained. Considering all the phases, a total of 6 equations 
for 6 unknown variables can be written. The resolution of this non¬ 
linear parametric system allows the magnitude and frequency of 
operation to be found. 

If we assume three identical phases, the condition of operation 
can be simply expressed as: 

- A three-phase double acting Stirling operation presents a 
120° phase angle between the motions of each piston and 
displacer. 

- The magnitude of the displacement of the piston of the 
engine i equals the magnitude of the piston of the engine 
i- 1 . 

Therefore, the following relationship is set between the input 
and output volumetric fluid flows Uj_i and u,- respectively: 

Ui = exp (-j»u,-i (35) 

In which tp = 120° for of a three-phase engine operation. 
Replacing Eq. (35), into Eqs. (33) and (34) allows the expression 
of the output pressure pf to be expressed as a function of the input 
pressure pf, such as pf = H pi (jVu)p) ,. The operation conditions are 
then: 

H P i(j(o)Z = - ^ (36) 

\H P ,(j(0)\ = 1 ( 37 ) 

The phase condition Eq. (36) determines the oscillation fre¬ 
quency co whereas the stroke is defined by the magnitude condi¬ 
tion Eq. (37). 

For T h = 120 °C, T c = 40 °C, the magnitude of the transfer func¬ 
tion H pi (jco ) at the 120° phase angle, is 1.027 at 35.36 Hz for 
0 mm stroke (case A, dotted line in Fig. 8). A limit cycle is not 
reached as condition Eq. (37) is not fulfilled. Still, the magnitude 
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Fig. 7. Simplified global network model of the engine i. 



Fig. 8. Magnitude and phase plots for for different strokes. 


of H pi (jco ) is greater than unity which indicates that the amplitude 
will grow with time. 

Higher strokes imply an increased stiffness of the beam at 
maximal amplitude and induce higher operating frequencies 
eventually. As an example, for a 1.5 mm stroke the frequency 
for which H pi (ja))Z = ^ is 44.33 Hz (case B, dashed line in 

Fig. 8). 

The limit cycle occurs at 57.35 Hz when the conditions of Eqs. 
(36) and (37) are reached simultaneously (case C, plain line in 
Fig. 8). In this case, the stroke is 2.08 mm. 

For higher stroke, the magnitude of the transfer function de¬ 
creases. For 2.2 mm stroke, the magnitude is 0.994 at 60.91 Hz 
(case D, dashed-dotted line in Fig. 8). This proves that a limit cycle 
is reached for the case C. 


3.2. Results and discussion 

3.2.1. Energy and power balance analysis 

Based on the electrical circuit of Fig. 7, the evaluation of the var¬ 
ious power quantities within the engine can be performed. Starting 
from the left hand side of Fig. 7, the average expansion space 
power (AESP) for one cycle is: 

it, = - / PdV e = 2 ReW i u ii) = u*) (38) 

In which u* , is the complex conjugate of the volumetric flow 
rate u,_! and Re is the real part. 

The same procedure can be used for the average compression 
space power (ACSP): 



AESP 2238.4 mW 


Phase 2 


Fig. 9. Power flow diagram for the 3 phases Stirling engine. 
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P 0 , = - jpdV c = ^Re(p c i ul i ) 


(39) 


For the predicted limit cycle (U p = 2.08mm at 57.35 Hz), the 
global power flow balance for the three-phase architecture (see 
Fig. la) is given in Fig. 9. The phases are numbered 1, 2 and 3. 

The power flow diagram of Fig. 9 describes the input and output 
power within the engine. As the engines are connected to each 
other, the AESP of the phase 3 equals the AESP of the phase 1. The 
total heat input can be seen as the AESP plus the regenerator losses 
which have to be compensated. The regenerator input power is 
equal to the total power generated by the engine (piston 
power + hysteresis loss + regenerator losses) as can be seen in Fig. 9. 

Therefore, the global efficiency is evaluated as: 

„ _ 281 + 161.2+16.7 _ -If) qccw 
I ~ 2238 . 4 + 16.7 — 

This value is the Carnot efficiency /y Comot = = 

20.35% which is in accordance with the isothermal assumption 
used in the model. 

The dissipated (useful) power at the piston of the engine 1 is 
281 mW which is 12.46% of the input power. 

The total losses sum up to 177.9 mW (7.89% of the input power) 
mostly associated to the gas spring hysteresis losses within the 
compression and expansion chambers (Eq. (31)). The ambient 
pressure operation, the low temperature conditions as well as 
the modest swept volume (31 cm 3 for a total volume of 440 cm 3 ) 
are in accordance with the small output power. 


3.3. Thermodynamics theoretical validation 


In order to perform a first validation of the proposed model, the 
Schmidt analysis is performed. This method is used in numerous 
studies [4,7] and its main results are recalled hereafter. From the 
ideal isothermal model assuming perfect gas law and given sinu¬ 
soidal motions of the piston and the displacer, the power output 
can be expressed as in Eq. (40): 


1 1 V swe k(t- l)sin(a) r ^ 

271 2fS Vsw yjx 2 + k 2 +2KTCOs(a) * ) 

(40) 

In which t, k and ft are defined in Eq. (41 ): 

T = — 

6 T h 

K = l “ ( 41 ) 

n _ \/2TKCOS(a)+T 2 +K 2 
P 2V+K+T 

Table 2 presents the comparisons between the Schmidt analysis 
and the network model for the same parameters used in Subsec¬ 
tion 3.2 (frequency and stroke U p ). No losses are considered in 
the Schmidt analysis. Therefore, the power evaluation used for 
comparison is the mean cycle power for one period. The results gi¬ 
ven in Table 2 clearly show that the two models matches. 


Pm = (OP,, 


3.4. Effect of the load 

The engine operating conditions can be driven tuning the load 
applied to the engine. In the present model, a single damping 


Table 2 

Comparisons between the linear network model and the equivalent Schmidt analysis. 



Schmidt analysis 

Network 

model 

Discrepancy 

(%) 

Power (mW) 

457.5 

459 

0.33 

Heat input (mW) 

2 247.4 

2255.1 

0.34 

Frequency (Hz) 

57.35 (given variable) 

57.35 


Stroke (mm) 

2.08 (given variable) 

2.08 



coefficient is used to model the dissipation and the load applied 
to the engine. We study here the impact of the two simultaneously. 

For T h = 140 °C, T c = 40 °C, the effects load on the stroke and the 
frequency are studied for a damping coefficient varying between 
0.4 N nr 1 s and 1.3 N s. Fig. 10a and b present the evolution 
of the stroke and the frequency respectively. 

From Fig. 10, it can be seen that low load cases lead to greater 
amplitudes and higher operating frequencies whereas the engine 
is slowed down with more moderate amplitude as the load is 
increased. 

As expected, the output power is enhanced for low load which 
is clearly shown by Fig. 11a. However, in this case, the efficiency is 
reduced because of the increase of the losses. Fig. lib plots the effi¬ 
ciency as a function of the damping. The efficiency is defined here 
as the power at the piston divided by the total heat input. Because 
the pressure drop losses lessen as the stroke decreased this effi¬ 
ciency increases with the damping. 

These results have to be read carefully. Indeed, the tempera¬ 
tures of the heater and cooler are supposed to be constant for 
the whole operation range. Still, these temperatures are strongly 
dependant on performances of the heat exchangers. 

Fig. 12 presents the evaluation of the Reynolds number for the 
fluid flow in the cooler, the heater and the regenerator. First, it can 
be noted that the results are in accordance with the proposed val¬ 
ues of the Reynolds number in Fig. 3. Second, the variations of the 
Reynolds numbers for the heater and the cooler are significant and 
the associated heat exchange coefficients are likely to be modified 
accordingly. Consequently, for given external temperature condi¬ 
tions, the internal temperatures of the engine will change and so 
will do the operating variables. 

The evolution of the temperatures as a function of the operating 
conditions (amplitude, frequency) and the geometrical parameters 
of the heat exchangers has been proposed in an analytical Stirling 
model in [21], For a given predicted operation point, such a model 
could be used to refine the results. Therefore, an iterative modeling 
strategy could be performed. The proposed network model to ac¬ 
count for the dynamic and the analytical Stirling model to take into 
account the heat exchanger performances can be used alterna¬ 
tively to converge to more accurate performance estimations. 


4. Experimental analysis 

A double acting three-phase Stirling engine has been built and 
is presented in Fig. 13. Its design is consistent with the details of 
Fig. 1. No pressurization is used and air is selected as the working 
gas. The heat input and cooling are not insulated to facilitate the 
instrumentation for the testing runs. 


4.1. Experimental engine description 

The heat source is obtained using 3 x 39 electric resistances 
connected in parallel. Therefore, 150 Watts heat power can be 
reached for each engine. However, in this first experimental test¬ 
ing, as no thermal isolation is used the heat power available for 
the engine is less than this value. 

For one engine, two temperature sensors are fixed on the back¬ 
sides of the cold and hot parts ( T C oLD_ex P and T H0T _exp respectively) 
whereas two other sensors are located within the compression 
and the expansion chambers (T c exp and T e _ exp )- Four additional 
temperature sensors are used to check the external temperatures 
of the two other engines. A pressure sensor is used and 3 acceler¬ 
ometers are fixed on the rigid links (see Fig. 13b). A constant elec¬ 
tric power is used during the experiment whereas the liquid 
cooling flow rate is fixed. 
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Fig. 10. (a) Evolution of the stroke with respect to the load, and (b) evolution of the frequency with respect to the load. 


(a) 




o 

cu 




Fig. 11. (a) Output power and losses for various loads and (b) relative efficiency-power parametric curve. 



frequency [Hz] 

Fig. 12. Evaluation of the Reynolds numbers for the heat exchangers and the 
regenerator. 


4.2. Experimental results and model comparison 

Fig. 14a presents the experimental stroke of the phase 1 at 
TcoLD_exp = 26.4 °C, T c exp = 22.1 °C, T eexp = 154.5 °C and T H or_exp = 


161.1 °C. The experimental operating frequency is 36.62 Hz 
whereas the stroke of the piston is evaluated at 0.15 mm. 

The assumptions used for the model are listed hereafter: 

- the temperatures of the model are set such as T e = T eexp and 
T c = T c _e X p ; 

- we assume that the three phases geometry and external con¬ 
ditions are the same; 

- the chambers of each engine are supposed to be symmetrical 
thus the volumes of the hot and cold chambers are equals; 

- the pressure drop is fixed at four times its theoretical value to 
account for the misalignment of the regenerator wire mesh; 

- the mass of the mobile component is estimated at 260 g; 

- the damping effect is only related to the membrane and used 
as the tuning variable. 

4.2.1. Steady state behavior 

The damping coefficient D of the model is modified to reach the 
operating conditions (Eqs. (36) and (37)). The value is set to 



Fig. 13. Experimental prototype, (a) Global view and (b) instrumented engine. 
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time [ms] [cm 3 ] 

Fig. 14. (a) Amplitude of motion of the piston membrane of the engine 1 and (b) Clapeyron diagram for the expansion chamber of the engine 1 for T cexp = 22.1 °C and 
T ejxp = 154.5 °C. 


4.33 Nrrr 1 s in which case the frequency is 36.42 Hz (0.56% differ¬ 
ence with experimental result) and the stroke is 0.15 mm (1.74% 
difference with experimental result). 

The experimental power at the expansion space (Fig. 14b) is 
18.3 mW whereas for the model it is evaluated at 11.26 mW which 
is 38.6% difference. The kinematic behavior matches though the 
model estimation of the pressure amplitude is not so accurate. 

These differences are thought to be associated to the symmetry 
assumption. Indeed, identified sealing flaws induced significant 
differences for the dead volumes of the three phases related to 
the different neutral positions of the membranes. Moreover, 
dissimilarities in the external temperatures have also been ob¬ 
served: phase 1 T C oLD_exp = 26.4 °C, T H 0 T_ex P = 161.1 °C; phase 2 
TcoLD_exp = 34.1 °C, T H oT_exp = 165.0 °C; phase 3 lcom_exp = 31.2 J C, 
THOT_exp = 163.6 °C. 

In such a case, the expansion and compression swept volumes 
values are not the same for a given phase contrary to the model 
assumption. The greater measured pressure amplitude in Fig. 14b 
may be produced by a compression swept volume greater than 
the expansion swept volume. 

4.2.2. Start-up conditions and self-starting property 

The experimental start-up temperature is evaluated at 
T c _ex P = 22.1 °C, T e exp = 145.1 °C whereas, for the model at the 
previously determined damping coefficient is 154.4 °C. Modifying 
the value of the linear damping to 4.092 Nnr 1 s allows the 
start-up temperature to be exactly estimated by the model. The 
difference between the damping coefficients is only 7.6%. 

The experimental self-starting of the engine whenever the 
critical temperature is reached is presented in Fig. 1 5. The engine 
is first stopped by clamping the membrane of the engine 1. Then 
the membrane is released and no external excitation added. A 



time [s] 


Fig. 15. Self-starting for T cexp = 21.7 °C and T eexp = 156.7 °C. 



4.37 4.38 4.39 4.40 4.41 4.42 4.43 


D [N m“'s] 

Fig. 16. Evolution of the stroke with respect to the damping coefficient. 


motion occurs within a short delay and the amplitude rapidly in¬ 
creases to a steady state value. The final stroke is 0.26 mm for a fre¬ 
quency of 36.07 Hz. 

It is worthy of note that in this case using the identified linear 
damping coefficient, the model gives a stroke of 0.66 mm for a fre¬ 
quency of 38.71 Hz. Increasing the damping coefficient allows to 
match the right amplitude. The required augmentation is only 2%. 

In order to obtain a first assessment of the stroke sensitivity to 
the damping in the experimental conditions, Fig. 16 shows the evo¬ 
lution of the amplitude with respect to the damping coefficient for 
the experimental proposed parameters and T cexp = 22.1 °C, 
T e _exp = 154.5 °C. It shows that for small strokes, the effect of the 
damping is the most important whereas it lessens for larger ampli¬ 
tude ( i.e. smaller damping). This high sensitivity is consistent with 
the identification of the different damping coefficient and its dras¬ 
tic effect on the behavior of the engine. 


4.2.3. Suggested strategies for optimization 

The engine has then been modified to improve the sealing and 
reduce the dissymmetry between the phases: 6 new silicone mem¬ 
branes have been realized especially. Consequently, the moving 
mass has increased to 310 g so the thickness of the beam springs 
has been modified to keep an operating frequency of about 
40 Hz. The new results are presented in Fig. 1 7. The experimental 
operating frequency is 39.98 Hz and the stroke adds up to 
0.426 mm for a slightly higher temperature difference. 

The new temperatures are set and the model is tuned adjusting 
the damping coefficient to 4.88 NitT 1 s (note that the dead volume 
have also been modified in the model to account for the new neu¬ 
tral positions of the membranes). Therefore, the theoretical fre¬ 
quency is 39.79 Hz (0.48% difference with experimental result) 
and the stroke is 0.424 mm (0.89% difference). The experimental 
power at the expansion space is 147.2 mW compared to the 
theoretical value of 108 mW. The power difference between the 
















































































F. Formosa et al./Energy Conversion and Management 78 (2014) 753-764 


763 




Fig. 17. (a) Amplitude of motion of the piston membrane of the engine 1 after modifications and (b) Clapeyron diagram for the expansion chamber of the engine 1 for 
T c _exp = 40.6 °C and T e exp = 178.2 °C. 


experimental and theoretical results is reduced from 38.6% to 
26.6%. It can also be noted that the output power has been drasti¬ 
cally increased. 

The improvement of the model prediction has been obtained 
while keeping the simplifying assumptions (symmetry of the 
engine, pressure drops in the heat exchangers neglected). With 
further modifications of the prototype to aim at a uniform heat dis¬ 
tribution between the phases the model-experimental correlation 
is expected to be improved. 

In the background of waste heat recovery, the performances of 
the engine would lie in efficiency and power output but also in 
standardization and economic criteria [3], The proposed simple 
model will be used to define a preliminary design for a given appli¬ 
cation. Some of the considered evolutions to be achieved are pre¬ 
sented hereafter and their effects will be easily assessed using 
the developed tools. 

Further studies will aim at the characterization of the heat 
exchangers and their optimization. To account for the exchanger 
and the regenerator heat transfer performances, the network 
model (which integrates the friction losses) will be exploited in 
an iterative strategy in conjunction with a previously developed 
thermodynamics model [21], The attractive solution which con¬ 
sists in the similarity of the hot and cold exchangers would also 
be assessed with respect to the performances gain with respect 
to standardization and cost reduction advantages. 

In order to operate the engine at higher pressure, the mechan¬ 
ical connections will be modified. Indeed, the planar proposed 
architecture would favor the interface with the heat source but 
lead to a cantilevered arrangement of the mechanical connections. 
The present engine cannot stand pressures more than one bar 
without detrimental large deformations. 

Additionally, the high nonlinear stiffness associated to the cho¬ 
sen beam structure would be lowered. The silicone membrane dis¬ 
sipation will also have to be reduced. By doing this, greater stroke 
amplitudes would be reached and larger swept volumes obtained 
eventually. The power of the engine will be enhanced accordingly 
providing that the pressure drop loss would be kept small. 

5. Conclusion 

The development of an equivalent electrical network model of a 
double acting FPSE has been detailed. The model allows the 
start-up condition to be established and the critical temperature 
evaluated. For given chambers temperatures and mean gas 
pressure, the steady state dynamic characteristics are estimated 
and the performances derived. A three-phase double acting free 
membrane engine architecture was proposed and is used as a test 
case. The amplitude of the engine is partly related to the non-linear 


stiffening behavior of the beam spring suspension which has been 
modeled and introduced in the model. An experimental engine has 
been built. The tests confirm the self-starting capability of the 
engine. The experimental results were compared to the model 
which is representative enough for further design works. The 
optimization of the heat exchangers will aim at improving the 
performances. Additionally, the mechanical connections and 
silicone membranes will be modified to allow for higher pressure 
operation and minimize the mechanical dissipation and the stroke 
amplitude limitation. From the ensuing results, the economic via¬ 
bility of double acting FPSE for waste heat recovery applications 
will be appraised. 
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